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The dynamics of a quantum system driven by an external field is well described by a unitary 
£T) | transformation generated by a time dependent Hamiltonian. The inverse problem of finding the 

field that generates a specific unitary transformation is the subject of study. The unitary trans- 
formation which can represent an algorithm in a quantum computation is imposed on a subset 



of quantum states embedded in a larger Hilbert space. Optimal control theory (OCT) is used to 

o 

CO 

o 
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state of Na2- Raman-like transitions through the A 1 ^^ electronic state induce the transitions. 
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solve the inversion problem irrespective of the initial input state. A unified formalism, based on the 
Krotov method is developed leading to a new scheme. The schemes are compared for the inversion 



Light fields are found that are able to implement the Fourier transform within a picosecond time 
scale. Such fields can be obtained by pulse-shaping techniques of a femtosecond pulse. Out of 
the schemes studied the square modulus scheme converges fastest. A study of the implementation 
of the Q qubit Fourier transform in the Na2 molecule was carried out for up to 5 qubits. The 
classical computation effort required to obtain the algorithm with a given fidelity is estimated to 
scale exponentially with the number of levels. The observed moderate scaling of the pulse intensity 
with the number of qubits in the transformation is rationalized. 



PACS number(s): 82.53.Kp 03.67.Lx 33.90.+h 32.80.Qk 
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I. INTRODUCTION 



Coherent control was initiated to steer a quantum system to a final objective via an ex- 
ternal field . If the initial and final objective states are pure, the method can be termed 
state-to-state coherent control. By generalizing, the problem of steering simultaneously a 
set of initial pure states to a set of final states can be formulated. Such a possibility has 
direct applicability in quantum computing where an algorithm implemented as a unitary 
transformation operating on a set of states has to be carried out irrespective of the input. 
In this application both input and output are encoded as a superposition of these states 

To implement such a control, the external driving field that induces a pre-specified unitary 
transformation has to be found. Different methods have been suggested for this task. Some 
rely on factorizing the algorithm encoded as a unitary transformation, to a set of elementary 
gates and then finding a control solution for the elementary unitary evolution of a single 
gate 0, 0] . The inherent difficulty in such an approach is that in general the field addresses 
many levels simultaneously. Therefore, when a particular single gate operation is carried 
out other levels are affected. This means that the ideal single gate unitary transformation 
has to be implemented so that all other possible transitions are avoided. The problem is 
simpler when each allowed transition is selectively addressable Q]. However, in general the 
problem of undesired coupling has to be corrected. A specific solution has been suggested 
(| but a general solution is not known. 

The presence of a large number of levels coupled to the external driving field is esp_e 
daily relevant in the implementation of quantum computing in molecular systems |7 
The use of optimal control theory (OCT) has been proposed as a possible solution 
In recent years OCT for quantum systems Q has received considerable attention, leading 
to effective methods for obtaining the driving field which will induce a desired transition 
between preselected initial and final states. To address the control problem of inducing a 
particular unitary transformation the state-to-state OCT has to be augmented. For exam- 
ple, if the unitary transformation is to relate the initial states \<fiik) with the final states 
\<Pfk), the state-to-state OCT derives an optimal field for each pair of initial and final 
states ((fik,<Pfk)- But the fields obtained are in general different, so that the evolution 
induced by e& is not appropriate for a different set of initial and final states. In order to 
implement a given unitary transformation a single field e that relates simultaneously to all 
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the relevant pairs (ipik,<Pfk) is needed. 

Two approaches have been suggested to generalize OCT for unitary transformations. The 
first approach is formulated directly on the evolution operator . An alternative approach 



121 ] . The present 



uses the simultaneous optimization of several state-to-state transitions jlO . 
paper develops a comprehensive framework for constructing an OCT solution for the uni- 
tary transformation. The study explores various approaches. A common framework for an 
iterative solution based on the Krotov approach 13] is developed. As a result, the numerical 
implementation of the methods are almost identical enabling an unbiased assessment. The 
implementation of the Fourier transform algorithm in a molecular environment is chosen 
as a case study. The performance of the various OCT schemes is compared in a realistic 
setup. A crucial demand in quantum computing is obtaining high fidelity of the solution. 
The present OCT scheme can be viewed as an iterative classical algorithm which finds a 
field that implements the quantum algorithm. The obvious questions are: 

• What are the computational resources required to obtain a high fidelity result? 

• How do these computational resources scale with the number of qubits in the quantum 
algorithm? 

• How do the actual physical resources i.e. the integrated power of the field scale with 
the number of qubits in the quantum algorithm? 

The paper is organized as follow: In Sec. |H] the problem is formulated introducing 
different objectives devoted to the optimization of a given unitary transformation. In Sec. 
Illll the application of the Krotov method of optimization of the given objectives is described. 
Expressions obtained for the optimal field are formulated as well as the implementation of the 
method. The variational method to derive the optimization equations is commented on in 
Appendix [XJ The results are used to study the implementation of a unitary transformation 
in a molecular model Sec. IIVI Finally, in Sec. |V| results are discussed. 
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II. IMPLEMENTATION OF A UNITARY TRANSFORMATION 



A. Description of the problem 

The objective of the study is to devise a method to find the driving field that executes 
a unitary transformation on a subsystem embedded in a larger Hilbert space. The unitary 
transformation is to be applied in a Hilbert space Ai of dimension M, expanded by an 
orthonormal basis of states (m = 1,...,M). The selected unitary transformation is 

imposed on the subspace J\f of N levels of the system (N < M) . In the context of quantum 
computation, the TV levels correspond to the physical implementation of the qubit(s) embed- 
ded in a larger system. The additional levels (m = N + 1, ...,M), considered as "spurious 
levels" , are only indirectly involved in the target unitary transformation. 

In any realistic implementation of quantum computing, "spurious levels" always exist. 
One reason is that the system is never completely isolated from the environment. In addition, 
the control levers, that in the present case is the dipole operator, connects directly only 
part of the primary levels. An example is the implementation of quantum computation 
using rovibronic molecular levels. The transition dipole connects two electronic surfaces 



[11 1 . The primary states reside on one surface, so that Raman-like transitions are used to 
implement the unitary transformation. The advantage of this setup is that the transitions 
frequencies between the electronic surfaces are in the visible region, for which the pulse 

n 

shaping technology is well developed 14]. Other levels residing on both of the electronic 
surfaces become spurious in the sense that any leakage to them destroys the desired final 
result. However at intermediate times these levels constitute a temporary storage space 
which facilitates the execution of the transformation. 

The objective is to implement a selected unitary transformation in the relevant subspace 
J\f at a given final time T. The target unitary transformation is represented by an operator 
in the Hilbert space of the primary system and is denoted by O. For N < M, the matrix 
representation of O in the basis has two blocks of dimension NxN and (M—N) x (M— 

N). The elements connecting these blocks are zero. This structure means that population 
is not transferred between the two subspaces at the target time T, but can take place at 
intermediate times. Only the NxN block is relevant for the optimization procedure, while 
the other remains arbitrary. 
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The dynamics of the system is generated by the Hamiltonian H, 

H(e) = H -/2e(t), (1) 

where H is the free Hamiltonian, e(i) is the driving field and /2 is a system operator 
describing the coupling between system and field. In the molecular systems, this coupling 
corresponds to the transition dipole operator and the driving field becomes radiation. In 
some cases more than one independent driving field can be considered. An example is when 
two components of the polarization of an electro-magnetic field are separately controlled 
The generalization of the formalism in such a case is straightforward. 
The system dynamics is fully specified by the evolution operator U(i, 0;e). An optimal 
field 6opt induces the target unitary transformation O at time T when 

U(T,0;e op ,) = e^O, (2) 

Eq. (J2J implies a condition on only the N x N block of the matrix representation of U. 
The phase (j)(T) is introduced to point out that the target unitary transformation O can be 
implemented only up to an arbitrary global phase. The phase can be decomposed into 
two terms, 0i(T) + <p2{T). The first, <f>x, originates from the arbitrary choice of the origin 
of the energy levels. Formally, a term proportional to the identity operator can always be 
added to the Hamiltonian. When the states \m) correspond to the eigenstates of H , the 



phase (pi is 



MT)= > (3) 



where E m is the energy of the level m. The phase 02 reflects the arbitrariness of the unitary 
transformation for the levels m — N + 1, M which are not part of the target. 

The method to determine the optimal field is based on maximizing a real functional of the 
field that depends on both the target unitary transformation and the evolution generated 
by the Hamiltonian, fulfilling Eq. (J2J). The problem of unitary transformation optimization 
is then reduced to a functional optimization. However different formulations of the problem 
can be made, leading to different functionals and then, in principle to different results. 
In the present context two different formulations have been proposed, one based on the 
evolution operator [ll] and the other on simultaneous iV state to state transitions jl^ . 
These formulations are closely related. A similar optimization procedure is described in Sec. 

nil 



B. Evolution operator formulation 



The optimization formulation is based on the definition of a complex functional r that 

n 

depends on the evolution operator at time T The following functional is introduced: 



N 

r(6; T; e) = Tr{6 f U(T, 0; e)P N } = ^(n|6 f U(T, 0; e)|n) 



(4) 



n=l 



where the projection P^v = Yun=i \ n )( n \ * s used. {|^)} denotes an orthonormal basis of the 
subspace Af. As O is a unitary transformation in the relevant subspace, the functional r 
is a complex number restricted to the interior of a circle in the complex plane of radius N 
centered at the origin. The modulus of r is equal to N only for an optimal field fulfilling 
Eq. (0). t can then be interpreted as an indicator of the fidelity of the implementation on 
the target unitary transformation 11]. When r approaches N, the transformation imposed 
by the field converges to the target objective. 

Since r is complex, several different real functionals can be associated with it. In Ref. 
[ll | the optimization of the real part of r, or the imaginary part, or a linear combination of 
both was suggested to find the optimal field. It was found that all these possibilities show 
a similar performance. For this reason, the present paper employs the optimization of the 
real part chosen as a representative case. The functional is therefore defined as: 



-Re 



r(0;T: el 



-Re 



N 



J>|6 T U(T,0;e)|n> 



n=l 



(5) 



The functional reaches its minimum value, F re = —N, when the driving field induces the 
target unitary transformation but however with the additional condition that the phase term 
exp(— i(f)(T)) is equal to one. 

Other functionals based on r but without any condition on the phase can be defined. In 
this work the squared modulus of r with a negative sign is studied: 



N N 



k(0:T:e) 



^TH6 f U(T, 0; e)\n) (n'\V(T, 0; e) t 6|n') 



(6) 



n=l rt'=l 



with minimum value F sm = —N 2 for any field satisfying Eq. 



C. Formulation of the Simultaneous N state to state transitions 



This formulation is based on the simultaneous optimization of N transitions between a 



set of initial states \l) and the corresponding final states 0\l) (I = 1, . . . , N) 



121. For this 



purpose the following functional is denned, 

{N ~\ N 

6 ] U(T, 0; 6) |0 (Z|U(T, 0; e)t6|i> = £ | (/|6 f U(T, 0; e) |/) | 2 . (7) 
i=i ) i=i 

Notice that while r is defined as the sum of amplitudes, 77 is defined as the sum of overlaps 

at the final time T. The parameter 77 is a positive real number and its maximum value, 

N, is reached when all the initial states \l) are driven by the field to the final target states 

0|Z), except for a possible arbitrary phase associated with each transition. The arbitrariness 

of these phases implies that the set of initial states \l) must be chosen carefully. In order 

to account for all the possible transitions, the states |Z) have to faithfully represent all the 

relevant subspace TV i.e. constitute a complete basis set. However, a choice of an orthonormal 

basis could produce undesired results. For example, the ambiguity of using an orthonormal 

basis {\n)} in the relevant subspace and an arbitrary unitary transformation D, diagonal in 

that basis. The product O D will also be a unitary transformation. If eo and €od are fields 

that generate O and O D at time T respectively, they both will have the same fidelity 77, 

r l± (d;T;eo) = V±(d;T;eo D ), (8) 

where _L denotes that 77 was evaluated using an orthonormal basis. Then any algorithm 
based on optimizing 77 that uses an orthonormal basis could find a solution for the field cor- 
responding to the implementation of an arbitrary unitary transformations of the form O D. 
(O is a particular case when D is the identity operator). The reason for this discrepancy is 
that 77 is only sensitive to the overlap of each pair of initial |/) and final 0\l) states, leaving 
undetermined the relative phases between them. For the optimization procedure to succeed 
a careful choice of the initial set of states is necessary. A simple solution is to compose the 
last N state as a superposition of all states in the basis Y2i=i \n)/VN, and to keep as is the 
first N — 1 states of an orthonormal basis. For this set of states, the maximum condition 
is achieved only when the field induces the target unitary transformation up to a possible 
global phase. 

To summarize the functional —77 in is used, 

F ss = - V (0;T;e), (9) 

with a minimum value F ss = —N. The optimal field reached satisfies Eq. (J2J), subject to a 
choice of the set of states |/) which determines the relative phases. 



D. Initial to final state optimization 



The present formulations of quantum control assume that the target unitary transfor- 
mation O is explicitly known, at least in the subspace N '. In most previous applications of 
optimal control theory, the objective was specified as the maximization of the expectation 



value of a given observable at time T subject to a predefined initial state Both mixed 
and pure initial states were considered jig. Il7|. A particular case is the determination of 
an optimal field to drive the system from a given pure initial state \tpi) to a target pure 
final state \(pf) at time T. This state-to-state objective optimization can be derived from 
the present formulation if the target unitary transformation becomes 0|</?j) = | </?/)■ The 
evolution operator formulation is then obtained by setting the projector = \<pi)(<pi\ in 
Eq. (HJ), obtaining the functional r, 

r(^, W T;e) = (<p f \U(T, 0; e)|^) . (10) 

The real functionals F re = — Re[r] and F sm = — |r| 2 are to be used in the study. As only a 
state-to-state transition is involved, the formulation is obtained by choosing \l) = \(pi). In 
this case i] = |r| 2 and F ss = F sm . Notice that this result is valid only when there is a single 
term in the sum in Eq. (@J) and in Eq. ((7j). 



III. OPTIMIZATION 



A common optimization procedure for all the functionals F as defined in the previous 
section is developed. The notation \n) for the states and n its index will be used in the 
evolution operator formulation. The notation \l) and I will be used in the simultaneous N 
state to state transitions formulation. The notation \tpik) and k where k = 1, . . . , N, will 
be used when the results are valid for both cases. An evaluation of any of the functionals 
requires the knowledge of the states \(pk(T)) = U(T, 0; e)\<fik) and \<ffk) = 0\<pik)- The 
operation of the evolution equation U(t, 0; e)\(fiik) can be calculated by solving the time- 
dependent Schrodinger equation 

j f \Mt)) = -^H(e)M*)>, (11) 

with an initial condition |y?fc(0)) = \<fik)- Since H = H(e) the state evolution will depend 
on the particular field. An alternative to Eq. (jllj) is the evolution equation for the unitary 
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transformation itself [ill ]. 

The method of optimization depends on the availability of the states of the system \<fk{t)) 
at intermediate times < t < T. 

Experimental realizations of OCT are typical examples where only initial and final knowl- 
edge of the states exist. For such cases feedback control and evolutionary methods are 
effective Such methods however require a large number of iterations to achieve conver- 
gence. A simulation of such a process requires repeated propagation of the N states by the 
Schrodinger equation. Thus they are computationally intensive. 

Computationally more effective methods are based on the knowledge of the states \<Pk{t)) 
at intermediate times. Additional constrains on the evolution are included that allow a mod- 



ification of the field at intermediate times consistent with the improvement of the o 



jjective 



18, 



the c_ .adient search m ethod Q, the Krotov m ethod Q, and the va— 
approach 21, A review of these common methods can be found in Ref. jl|. In the 
present study the Krotov method has been adopted. A brief description of the alternative 
variational method is given in Appendix |XJ 



A. Krotov method of optimization 

The Krotov method is utilized to derive an iterative algorithm to minimize a given func- 
tional that depends on both final and intermediate times Cf. Ref. j^. 

For convenience, the equations are stated using real functions: atfc m (t) = Re[(m|^(*)}] 
and Pkmif) = lm[(m\ipk(t))}. The notation cx k and f3 k is used to describe the M-dimensional 
vectors with components a km and (3 km . Using such a notation, the evolution equation (fTTj) 
becomes: 

-cc k {t) = Sl R (t, e) • a k {t) - e) - (3 k (t) , 
at 

j t f3 k (t) = nj(t, e) ■ a k (t) + n R (t, e) ■ (3 k (t) , (12) 

where and Qi are real matrices with the corresponding components composed of the 
real and imaginary parts of Clij = (i\(—iH/h)\j} where \i) and \j) are states from the basis 
set The initial conditions are given by the vectors a. k (0) and (3 k (0) with components 
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composed of the real and imaginary part of the amplitudes (m\ipik)- cxjk and f3j- k denote 
the vectors corresponding to the amplitudes (m| </?/&). 

The formalism considers a. k , f3 k , and the field e to be independent variables. A necessary 
consistency between them will be required in the final step of the algorithm. The vectors 
f a and fa constitute the right hand side of Eq. 1)1 2|). 

f a (t,a k ,P k ,e) = ft R (t,e) ■ a k (t) - Sl^t, e) ■ (3 h {t) , 

f p {t, a h , (3 k , e) = U T (t, e) ■ a k (t) + Q R (t, e) • (3 k (t) . (13) 

The vectors f a (fa) are equal to the total time derivative of a (f3) only when the state 
is consistent with the field through the evolution equation (JT2J). The dependence of a and 
(3 on t will be made explicit only when necessary. An important property of the problems 
under study is that f a and fa are linear in the functions {a, /3} and e. This choice simplifies 
the optimization problem, the non-linear case has been studied in Ref. 23^ . 

A "process" w = w[t, {et, f3}, e] is defined as the set {a, f3} of iV vectors cx k and iV vectors 
/3 k related to the field e through the evolution equations with the initial conditions CKfc(O) 
and /3fc(0). A functional of the process can be defined: 

J[w] = F({cx(T),(3(T)}) + f g(e)dt. (14) 

Jo 

For the present applications F can be any of the functionals F re , F sm , F ss as introduced 
in the section |H] The optimal field is found by a minimization of the functional J. The 
integral term represents additional constrains originating from the evolution equation of the 
system. For simplicity only the case where g is a function of the field e is presented, but a 
generalization to a more general case in which g depends on ot(t) and f3(t) is straightforward. 
The particular dependence of g on the field will be discussed later. 

The main idea in the Krotov method is to introduce a new functional that mixes the 
separated dependence on intermediate and final times in the original functional (|14j). Using 
the new functional it is possible to derive an iterative procedure that modifies the field at 
intermediate times in a consistent way with the minimization of F at time T. The new 
functional is defined as 

L[w,$\ = G({ct(T),(3(T)}) - $(0,{«(0),/3(0)}) - f R(t,{a, /3},e) dt , (15) 

Jo 

where, 

G({ a (T),(3(T)}) = F({a(T), (3(T)}) + $(T, {a(T), /3(T)}) , (16) 
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and 



<9$ 

R(t,{a,P},e) = -g(e) + — (t,{a,/3}) 
N d<f> 



+ Y] q — [t,{a,/3}) ■ f a (t,a k ,p k ,e) 



k=l 

{a,/3}) denotes an arbitrary continuously different iable function. The partial deriva- 
tives of $, and J|J^, form a vector with m components. In the following t, a and /3 are 
considered to be independent variables in $. 

When {a, f3} and the field are related by Eq. (|12|), R can be written as R = —g + d<&/dt. 
Introducing this result in Eq. (|15|). it can be shown |23[ that for any scalar function $ and 
any process w, L[w, $] = J[w]. Then the minimization of J is completely equivalent to the 
minimization of L. 



1. Iterative algorithm to minimize L 

The advantage of the definition of the functional L[w] is the complete freedom in the 
choice of $. This property is used to derive from an arbitrary process [t, {a.(°\ (3^}, e^} 
a new process w^[t,{a^,P^y,e^] such that L[w^\^] < $]. The procedure can 

be summarized as follows: 

• (i) $ is constructed so that the functional L[«/°)] is a maximum with respect to any 
possible choice of the set {a, f3}. This condition gives a complete freedom to change e. 
The related changes of the states are consistent with the system evolution, Eq. (|T2|). 
and therefore, will not interfere with the the minimization of L. 

• (ii) A new field is derived with the condition of maximizing R, decreasing then the 
value of L with respect to the process w^°\ In this step the consistency between the 
new field and the new states of the system {a^\ p* 1 '} must be maintained. 

The new field becomes the starting point of a new iteration, and steps (i) and (ii) are 
repeated until the desired convergence is achieved. 
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2. The linear problem: construction of $ to first order 



The difficult task in the Krotov method is the construction of $ so that L is maximum 
for {a(°\/3 (0) }. The maximum condition on L is equivalent to imposing a maximum on G 
and a minimum on R. However, in some cases the maximum and minimum conditions can 
be relaxed to extreme conditions for G and R, which simplifies the problem. 

The extreme conditions for R with respect to {et(°\(3^} are given by 

9fi (i> ( ^ (0) } )£ (0) ) = 0, 

(t,{a(°),^ )}, e (°)) = 0. (18) 



dct k 
dR 



d/3 k 

The following vectors are introduced 

7*(t) = !^(U« (0) ,/3 {0) }), 

OOL k 

Suit) = ^-(t,{^°\(3^}). (19) 
o/3 k 

y k and S k are only functions of t, as the partial derivatives are evaluated in the specific set 
{a.(°\ f3^}. Using Eq. (|T2j) the extreme conditions can be written as 

j t 6 k (t) = nftt, e W) • lk {t) - n T R (t, e (°)) • s k (t) , (20) 

where f2 T denotes the transpose of the matrix f2. The extreme conditions for G are 



9a fe (T) 



({«W(T),/3^(T)}) = 0, 



/9G 

— -({Q(°)(T),^(r)}) = 0. (21) 
oph\T) 

Using Eq. JEIl and Eq. JEJ) 

7*C0 = -^-7^({«(0)(T),/3 (0) (T)}), 

$k(T) = -— — ({a<°>(T),/3(°>(T)}). (22) 
dp k (T) 

The above conditions at time T, together with Eq. (|2U|) determine completely the set 
{jit), S(t)}. As they are defined as the partial derivative of $ with respect to a km and (3 km , 
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$ is expanded to first order (denoted as $*), 

N 

$*(*, {a, /3}) = £ [ 7fc (0 • a*(f) + <5 fc (t) • fc (f)] . (23) 
fc=i 

By employing $*, the functions G* and i?* can also be constructed to first order using Eq. 
()16|) and ()17|) respectively. This completes the first step in the iterative algorithm. 

To accomplish the second step R is maximized with respect to the field. Again in some 
cases the maximum condition can be relaxed to the extreme condition dR/de = 0. Using 
the expression for R* leads to 

k=l 

+ E«.(«)-^(*,^ w ,i8f,««). (24) 

k=l 

Eq. ()24|) is used to derive the new field e 1 ) in each iteration. This equation must be solved 
in a consistent way with Eq. (jl2J) describing the system dynamics. 

Due to the use of extreme instead of maximum or minimum conditions, it must be checked 
that the new process, improves the original objective in each iteration J[w^] — J[w^)] > 
0: 

J[w (0) ] - J[w (1) ] = L[#,f*] -L[wM,$*] = At + [ A 2 (t)dt. (25) 

Jo 

where, 

Ai = G*({ct {0 \T),f3 (0) (T)}) - G*({a«(T),/3 (1) (T)}). (26) 

A 2 (t) = I?b{cP\ - tf{t,{ a M,0P>},eV>). (27) 

The above relation is obtained when f a and fa are linear in {ot,[3} 

R*(t,{ a ,f3},e^) = -g(t,eM), (28) 

for any value of {a, f3}. 

A sufficient condition for Jfw*- -*] — Jf-u/ 1 **] > is Ai,A 2 (t) > 0. Ai depends on the 
choice F and A 2 (i) on the choice of g so that each case must be analyzed separately. These 
conditions imply that the Krotov iterative algorithm convergence monotonically to the final 
objective. 
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3. Dependence on F 



The dependence of G* on F can be made explicit by introducing $* and using Eq. 
and Eq. (flfijl . 

G*({cx(T),(3(T)}) = F({cx(T),(3(T)}) 

" OF dF 
^ da k (T) d(3 k {T) 

When F is linear in {a, (3} G* = and then Ai = 0. In this case all the improvement 
towards the original objective in the iteration is due to the term Aa(t). When F is non 
linear in {a, (3} the condition Ai < must be checked in each case. 

An additional difficulty is that the conditions (fT9*|l for 7 and 8 depend on the particular 
choice of F. In all the cases under study (F re , F sm and F ss ) 

7* CO = c k a fk {T), 

8 k {T) = d k (3 fk (T), (30) 

where the coefficients c k and d k depend on the sets {a^°\T), /3^°'(T)} and {otf,f3f}. Defin- 
ing the vectors 7 fc = c k ~ 1 j k and 8 k = d^ 1 8 k , the conditions (J22J) for all the cases under 
consideration are: 

%{T) = a fk {T), 

~6 k (T) = (3 fk (T). (31) 

Their evolution is given by Eq. (|2()|). Changing 7 and 8 to 7 and 8 Eq. (|24jl can be written 

as 

fc=i 

+ f:4^(t)-^(t,««,/3«, e W). (32 ) 

k=l 

The different choices of F imply different coefficients (c k and d k ) and a possible different set 
of initial \<Pik) and final \(pfk) states. Nevertheless, the iterative procedure is identical in all 
the cases. 
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4- Dependence on g(e) 



A delicate point is the choice of the function g(e) in J[w]. The time integral in the func- 
tional should be bounded from below, otherwise the the additional constraint will dominate 
over the original objective F in the functional J. In addition A 2 (£) > is required in order 
to guarantee the monotonic convergence of the optimization method. 

A consequence of the linear dependence of f a , f p and Eq. is that the function R* 
for the new process has the simple form 

R*(t,{ccM,(3M},eM) = -g(e) + (e« - e <P>) |(e«) . (33) 

Using this expression together with Eq. (J2H|) leads to 

A 2 (t) = -g(e^) + g(e^) + (e« - e<°>) |(e«) , (34) 

A choice of g(e) fulfilling the previous requirements is 

g(e) = X(t) [e(t) - e(t)] 2 , (35) 

where e is a reference field and A(t) is a positive function of t. Using Eq. (j3^|) and for any 
field e 

A 2 (t) = \{t) [Ae(t)] 2 > 0, (36) 

where Ae(t) = e^(t) — e^(t). The method therefore will converge monotonically. Using 
Eq. and Eq. (|33|) the field in the new iteration becomes: 

N 

= e(t) LVJ„ + ) ■ '-J-lli rv l[] r'H- 



+ d k ~S k (t).^l(t,a^,^\e^)\ . (37) 



de 

The result of the iterative algorithm depends strongly on the choices of the reference field e 
and on the function \(t). 

Two possible choices of e are analyzed. The first, e = 0, is the one commonly used 
in OCT applications Q]. In this case, the additional constraint in J[w] has the physical 
meaning that the total energy of the field in the time interval [0, T] is limited. This however 
presents a problem when the iterative procedure reaches the optimal field. The iterative 
method is found to reduce the total objective J by reducing the total pulse energy, slowing 

15 



and even spoiling the convergence to the original objective F. The usual remedy is to stop 
the iterative algorithm before this difficulty is reached. However, such a procedure could 
prevent the optimization algorithm from obtaining high fidelity. 
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16, 231. In this iterative 



A different possibility is e = e^ -* can avoid this problem 
algorithm must be interpreted as the field in the previous iteration. Now the additional 
constraint in J[w] has the physical interpretation that the change of the pulse energy in 
each iteration is limited. When the iterative procedure approaches the optimal solution 
the change in the field vanishes. Therefore, the convergence to the original objective is 
guaranteed. In the rest of the study e = was chosen. 

The function \(t) introduces the shape function s(t) i.e. \(t) = \ /s(t). The purpose 



m. A 



of s(t) is to turn the field on and off smoothly at the boundaries of the interval 
is a scaling parameter which determines the optimization strategy. When Ao is small the 
additional constraint on the field in the functional becomes insignificant, resulting in large 
modifications in the field in each iteration. This is equivalent to a bold search strategy where 
large excursions in the functional space of the field take place with each iteration. Large 
values of Ao imply small modifications in the field in each iteration, slowing the convergence 
process. Using large values of A is a conservative search strategy which is advantageous 
when a good initial guess field can be found. A possible mixed strategy is to initially use a 
bold optimization with small values of Ao- This leads to a guess field for a new optimization 
with a large value of An 



B. Application to the functionals F re , F sm and F s . 



Based on the derivation of the Krotov method it is possible to connect directly the 
minimization F re , F sm and F ss to the correction to the field. Eq. (|2T)|) corresponds to the 
evolution of a set of states {\Xk(t))}, 



d 



Xk(t)) 



H 



\Xk(t)) 



(38) 



dt h 

with the conditions ()31|) . \xk{T)) = \tpfk)- The formal solution of the equation is given by 
\Xk(t)} — U(t, T; e^)0\(fik)- Using Eq. (|3T|) the correction to the field in each iteration 
becomes: 

S « Im 

_fc=l 



Ae(t) 



An h 



N 



^(e {0) ) & ik \d j v\t, T; e (°>) Jl U(t, 0; e^)\cp ik ) 



(39) 
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The coefficients a&, will depend on the particular choice of the functional F and are related 
to Cfc and dk defined in Eq. (JT^J. For F re and F sm , the states {|y?ifc)} denote an orthonormal 
basis {|n)} of the relevant subspace Af. For F re the coefficients are a n)Te = 1/2, and as this 
functional is linear on the states Ai re = 0. For F sm the coefficients are 



N 



(ri\TJ(p,T;e<®)6\ri) 



(40) 



n'=l 



and then are equal for all the states in the basis of N '. In addition, 



Ai, 



N 

X>| (u(T,0;6(°))-U(T,0;e«))6|n) 

n=l 



(41) 



Therefore, A l sm > 0. For the F ss functional the set {\(Pik)} for which the states are denoted 
by {\l)} the coefficients a>i are 



ai,. 



<i|U(0,T;e<°>)6|Z> 



(42) 



depending on the index I corresponding to each state. In this case 



N 

A llM = (u(T,0;e(°))-U(T,0;eW))6|0 



(43) 



and then Ai jSS > 0. 

The results Ai > and ^(t) > guarantee the monotonic convergence of the iterative 
algorithm based on the Krotov method for the three functionals. 



C. The optimal field 



The optimal field has the property that the field correction in the next iteration Eq. (|39|) 
should vanish. Defining this correction as: 



C(t;e) = Im 



a ^ (0) ) Wo'u'ft T; e) £U(t, 0; e)|^) 



fc=i 



(44) 



where e is an arbitrary solution for which C(t;e) = 0. 

The first question to be addressed is whether any optimal field, defined by Eq. (J2J), is 
a possible solution of the iterative algorithm. e opt denotes a field that generates the target 
unitary transformation up to a global phase, U(T, 0; e opt ) = e~ 1 ^ O. Using the relation 



U(t,0;e)=U(t,T;e)U(T,0;e). 



(45) 
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In addition, the relation \^/ k (t)) = U(t, T; e opt ) 0\(pn-) implying that the term 
k(t)\fl\^ k (t)) is real, simplifies Eq. (jSJ) to: 

jv 

C{t-e opt ) = ^(# fc (t)|/2|* fc (t))Im 

fc=i 

Using Eq. (J4"UJ) for the functional F sm leads to a„ iSm = iVexp in Eq. (|44|). A similar 
result is found for the functional F ss , a; iSS = exp given by Eq. (|4"2j). Therefore any field 
generating the target unitary transformation is a possible solution of the iterative algorithm 
based on any of the functionals F sm and F ss . This result does not imply that when initializing 
the different iteration schemes with the same guess field the same solution will be obtained. 

The analysis is more complex for the functional F re . The coefficients are now 
real a n>re = 1/2, and are independent of the state index. This leads to C re = 
Im[exp (— i<f>)] ^2 n= i{& k{t)\ P- \^k(t))- The sum is generally different from zero and the so- 
lutions to the algorithm are fields with a phase term exp (— iifi) = ±1. Only the case +1 
minimize the original functional F re , but the relaxation to extreme conditions in the Krotov 
method allows to obtain other physically valid solutions. In the special case in which the 
unitary transformation is imposed on all the Hilbert space (N = M), any optimal field is a 
possible solution regardless of the global phase. The reason is that the sum in C re is zero 
since /z is a traceless operator. 

The phase sensitivity of the functional F re can be demonstrated in the state-to-state 
optimization. The iterative algorithm in this case will converge to a field that drives the 
system to the final state +\<pj) or —\(pf), while the optimization of F sm or F ss will converge 
to the final state up to an arbitrary global phase. There is no a priori advantage however 
to any of the three functionals in the convergence rate or in the simplicity of the solution. 
The solutions are physically equivalents since they differ only in a global phase. 

In addition to the desired optimal fields, the algorithm could also generate spurious 
solutions. A possible example is the functional F sm employed to implement a unitary trans- 
formation Od with a matrix representation diagonal in the basis of the free Hamiltonian 
eigenstates \e n ). In such a case C(t,e = 0) is proportional to the diagonal matrix elements 
(e n \/2\e n ). When these matrix elements are zero, e = is a solution of the iterative algo- 
rithm, but it does not implement the desired unitary transformation. A simple remedy to 
overcome this difficulty is to use a different initial guess to start the algorithm. 



(46) 
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D. Discrete implementation of the optimization algorithm 



A numerical solution of the iterative optimization algorithm requires a discretization 
scheme for the time axes. The correction to the field Ae is implicit and appears on both 
sides of Eq. (|39|) . To implement the procedure, two interleaved grid points in time were 
used. The first grid was used to propagate the states. The second grid was used to evaluate 
the field. The grid describing the states has N t + 1 points separated by At = T/N t , from 
t = to t — T. The grid representing the field has Nt points separated by At and starting 
at t = At/2. The initial set of states \(fik) was used for the target unitary transformation 
O optimization with the functionals F re} F sm or F ss . The numerical implementation of the 
algorithm follows: 

• (i) Using an initial guess field , the states (p fk are propagated in reverse from t = T 
to t — to determine U(t,T; e^)0\ipik) on the time grid of states. 

• (ii) The new field is determined in the interleaved grid point t = At/2 using the 
approximation 

N 



MA( /2)«-^to 

An h 



*k(e {0) ) (^|6 T U T (0, T; e<°>) fi U(0, 0; e (1) )|^ fe ) 

-fc=i 

(47) 

Notice that U(0, 0; e^ 1 )) |^.) = \(pik)- Then the new field in the first field time grid 
point is obtained, e <yl \At/2) = e' ' + Ae(At/2) and used to propagate \(fik{t = 0)) to 
the next state grid point t = At. The same process is used to obtain the new field 
e 1 ) in the next field time grid point t = At + At/2, evaluating the correction with the 
already know states in the state grid point t = At. The process is repeated to obtain 
e^ 1 -* in all the field time grid points. 

• (iii) The new field e*- 1 -* is used as input to the new iteration (e 1 - -* = e*- 1 -*) and the process 
is repeated until the required convergence is achieved. 

More elaborate methods to deal with the implicit time dependence of Eq. (|39|) have been 
developed. For example, approximating the dynamics in between grid points by the free 



evolution with H r 



221 ] . The simple procedure, which is able to keep the monotonic behavior 



of the optimization method was found sufficient. 
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The present implementation is based on a forward time propagation. Using the same 
formalism, the optimization can be accomplished also by a backward time propagation. It 
is also possible to combine both cases, and to perform the optimization in the forward and 
backward propagations 0, |3| . In the current studies, these other procedures were found 
to be inferior, slowing down the convergence rate. 



IV. THE FOURIER TRANSFORM EXAMPLE IN A MOLECULAR MODEL 

As an illustration the implementation of Q qubit Fourier transform in a two electronic 
surfaces molecular model was studied. Fig. ^ shows a schematic view of a model based on 
the electronic manifolds of Na2. 

The Hamiltonian of the system describes a ground and excited electronic potential energy 
surface coupled by a transition dipole operator: 

H = H 9 <g) \G)(G\ + H e ®\E)(E\ -ft® (\G)(E\ + \E) (G\) ■ e(t) (48) 

where \G) and \E) are the ground and excited electronic states and H 9 and H e are the cor- 
responding vibrational Hamiltonians. The electronic surfaces are coupled by the transition 
dipole operator £t, controlled by the shaped field e(t). 

The present model is a simplification of the Na2 Hilbert space restricting the number 
of vibrational levels. On the ground X electronic state the first 40 vibrational levels 
selected out from the 66 bound states are used. In the excited state, the lowest 

20 vibrational states are used out of the 210 bound levels. The vibrational Hamiltonians 
become therefore 

40 20 

H g = J2 E gi \gi){gi\ ; H e = J2 E ej\ e j)( e j\ ■ ( 49 ) 

For Na2 the 00 transition frequency between the ground vibrational levels of each surface is 
Q = E el — E g i « 0. 06601a. u. (~ 1.8eV ). A transition dipole operator independent of the 
internuclear distance R was considered, p, = fi (\G)(E\ + \E)(G\). This model is sufficient 
for the illustrative purpose of demonstrating the execution of an algorithm in a molecular 
setting. 

The N = 2® first levels of the ground electronic surface are chosen as the registers repre- 
senting the Q qubits. The unitary transformation implemented is a Fourier transform [^J 
invoked on the iV levels on the X l Y^ electronic state representing the qubit(s). The unitary 
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3 5 7 9 11 13 15 

R(a.u.) 

FIG. 1: Schematic representation of a molecular model based on the vibrational levels in the X^Y^ 
(lower) and -A 1 £+ (upper) electronic surfaces of the molecule Na2- Atomic units are chosen h = 1. 
R denotes the internuclear distance. The arrows indicate two of the possible transitions induced 
by the driving field between arbitrary levels in the lower and upper surfaces. On the right is a 
magnified view of some of the energy levels involved and transitions between them. 

transformation is implemented through transitions between the two electronic manifolds Cf. 
Fig. HI 

An implementation of the iterative algorithm is chosen where the |^) <g> \G) and \ej) <8> \E) 
eigenstates are used as the basis The N = 2® first states in the lower surface are 

used as the basis {\n)} of the relevant subspace. The first N — 1 energy levels plus the 
linear combination J2n=i \n)/VN are used as the set |Z) for the state to state formulation. 
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The wavefunction propagations were carried out by using a Newton polynomial integrator 
. The final time for the implementation is T = 4.5 x 10 4 a.u. (« 1 psec). In all the 
cases a Gaussian shape function s(t) = exp{— 32(t/T — 1/2)} and a guess field e guess (t) = 
e s(t) cos(fit) were chosen. 

0.00 
-0.25 
E -0.50 

o 
c 

~i 

-0.75 
-1.00 

50 100 150 

iteration 

FIG. 2: Normalized functional, J n orm versus the number of iteration: F re (squares), F sm (circles), 
F ss (triangles up) for implementing a FFT in 4 levels. The line with triangles pointing down 
corresponds to F ss functional when {\l)} is chosen as the orthogonal basis The objective is 

reached when J n0 rm = — 1- Ao = 10 3 and eo = 5 x 10~ 3 a.u. in all the cases. 

The implementation of the Fourier transform in 2 qubits (N = 4) embedded in the set 
of 60 levels is used for comparing the performance of the methods. Fig. El shows the change 
in the normalized functional, defined as J norm = J/N for F re and F ss , and J norm = J/N 2 
for F sm , with the progression of the iterative algorithm. In all the cases the target value 
of the normalized functional is —1. A large reduction in the value of the functionals is 
accomplished in a small number of iterations. Notice the behavior of the simultaneous state 
to state formulation F ss with an insufficient choice of the states The algorithm finds a 
minimum of the objective, but, as shown in Fig. 0J the fidelity saturate at a very low value 
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meaning that this field does not generate the target unitary transformation. 




FIG. 3: Evolution of the optimization in the complex r plane for the case in Fig. [2j The lines 
correspond to F sm (circles), F re (squares), and F ss (triangles up). The open circle indicates the 
value of r for the common guess field. The dashed black line is the circle |t| = N indicating the 
target of the methods. The arrows mark the direction of convergence. The insert enlarges the 
region corresponding to the real axes close to the circumference. 

Fig. El shows the value of r for the field obtained in each iteration. The same initial 
guess was used in all the cases which constituted the starting point for all the iterative 
optimizations. However, the final results depend on the particular functional used. As 
discussed before the method based on F re finds a solution with a phase factor exp (—i<p) ~ 
+1. 

For the purpose of quantum computing the target unitary transformation has to achieve 
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iteration 

FIG. 4: Fidelity of the implementation of the 2 qubit Fourier transform versus the number of 
iterations N%t for the optimization in Fig. [2j The lines correspond to F re (squares), F sm (circles), 
and F ss (triangles up). F pop (triangles down) denotes the case when the set {\l)} is chosen as the 
orthogonal basis {\n)} for the functional F ss . 

high accuracy. The fidelity functional 

fidelity = log 10 (l - \r\ 2 /N 2 ) . (50) 

is used to indicate the quality of the solution. Fig. 0] shows the improvement of the fidelity 
versus the iteration. The square modulus functional F sm Eq. © shows a faster convergence 
rate than the other two functionals. 

In Fig. 03 the Fourier transform of the field for each of the optimization procedures is 
shown. The large peak at the frequency Q, seen in all cases, is the result of the choice of the 
guess field. Besides a similar width in frequencies is found. However, the fidelity reached 
by the solution corresponding to the square modulus functional F sm is significantly better 
than in the other cases for the same number of iterations. 

The molecular model is also used to compare the convergence of the unitary transforma- 
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CO 

FIG. 5: Fourier transform of the optimal field result of the optimization in Fig. 0]for the functionals 
Frej ^imi and F ss . 

tion with the size N of the relevant subspace. Fig El shows the improvement in the fidelity 
versus the number of iterations for implementing a Fourier transform in 2, 4, 8, 16, and 32 
levels (1, 2, 3, 4, and 5 qubits respectively). The convergence characteristics in the initial 
iterations strongly depends on the initial guess and the parameter Ao- For example the 
initial guess seems inappropriate for the 1 qubit case which displays an initial very slow con- 
vergence until after 25 iterations it find the right track. After a large number of iterations 
the convergence characteristics settled meaning that each new iteration was only a slight im- 
provement on the previous one. As the iteration proceeds the rate of convergence decreases 
in all cases, scaling approximately as the inverse of the number of iterations. Comparing the 
rate of convergence for the different number of qubits after a large number of iterations the 
rate seems to be inversely proportional to the number of levels. High fidelity was obtained 
for 1, 2, 3 qubit cases by continuing to 600 iterations. The results allow to compare the 
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FIG. 6: Fidelity versus the number of iterations for implementing a Fourier transform in 2 (dashed- 
dotted line), 4 (circles), 8 (solid line), 16 (dashed line) and 32 (dotted line) levels. 

integrated intensity of the optimal field: 

1=1 W(t)\dt (51) 
Jo 

The initial integrated intensity for all cases was identical. The optimization procedure 
changed X depending on the number of qubits. The converged results show a moderate 
increase of X with the number of levels starting from X = 42 for 1 qubit to X = 54 for 2 
qubits and X = 78 for 3 qubits. 



V. DISCUSSION 



An implicit assumption in the optimization procedure is that the system is controllable. 
This means that a field e(t) exists which implements the unitary transformation up to 
a pre-specified tolerance. The problem of controllability has been the subject of several 
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32 



331 ] . In the context of unitary transformations it has been shown 



studies |2I 

(^(j] that if the commutators of the operators H and /2 generate the complete Lie group 
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SU(N), the system is completely controllable. In more concrete terms addressing the Na2 
model, it is expected to be completely controllable. The reason is that the energy levels 
are non degenerate and in addition each transition is distinct, characterized by a different 
Frank Condon factor (ej\fj,\gi). This controllability property will be true in almost any 
non-symmetric molecular system. 

A far reaching conclusion is therefore that for any unitary transformation contained in the 
Hilbert space of the molecule, there is a driving field that implements the transformation in 
one step. In a molecular system this task could be achieved in a time scale of a picosecond. 
Since a field that executes such a unitary transformation exists, how difficult is it to find it? 
Does this optimal field have reasonable intensity and bandwidth? 

The OCT scheme can be considered as a classical algorithm employed for the inverse 
problem of finding the field that generates a predefined unitary transformation. The diffi- 
culty of the inversion process is related to the scaling properties of the numerical effort with 
respect to the number of levels N. The best OCT algorithm based on the F sm functional is 
then used for estimating the scaling. 

Simulating the quantum evolution is the major numerical task of the algorithm imple- 
menting OCT. The basic step is a single vector matrix multiplication which represents the 
operation of the Hamiltonian on the wave function. This task scales as 0(M 2 ) for direct 
vector-matrix multiplication or 0(M log M) for grid methods based on FFT 3^|. The time 
propagation requires Nf steps which scale as 0(TAE), where AE is the energy range of the 
problem. 

The simulation of a unitary transformation in the relevant subspace turns out to be N 
times more costly. Summarizing, the numerical cost of the classical simulation of the quan- 
tum propagation scales as Cost ~ 0(2® M 2 TAE). This scaling relations is consistent with 
the fact that a classical simulation of a quantum unitary transformation scales exponentially 
with the number of qubits. 

The numerical cost of the OCT iterative algorithm used for inversion can now be ex- 
amined. The crucial question is how many iterations are required to obtain the field that 
implements the unitary transformation up to a specified fidelity /. The analysis of the re- 
sults of Sec. IIVI show that only the initial iterative steps are very sensitive to the choice of 
the initial guess field. Eventually the rate of convergence reaches an asymptotic behavior 
where the fidelity becomes inversely proportional to the number of the iterations steps. In 
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addition Cf. Fig. 13 the rate of convergence is also inversely proportional to the number of 
levels. This relation implies that the number of iterations Nu required to achieve the fidelity 
/ becomes 

N lt « be^ , (52) 

where the coefficients a and b are positives. The data confirm that the coefficient a is 
independent of the number of levels N. The consequence of Eq. (|52|) is that the numerical 
resources required on a classical computer in order to implement the proposed scheme, scale 
exponentially with the number of levels N. Finding the field that implements in a single step 
a large unitary transformation is therefore prohibitively expensive. Thus fields that achieve 
high fidelity are only feasible for unitary transformations with a small relevant subspace. 
The limiting case would be the one dimensional state-to-state optimization. 

Quantum control is based on interferences between many distinct pathways 1]. State to 
state coherent control finds a constructive interference which leads exclusively to the final 
state. The controllability depends on having a sufficient amount of interference pathways. 
Implementing a unitary transformation by interferences is more complex. In this case the 
interference pathways from one state to another have to avoid other interference paths which 
connect other states. The possible number of interference pathways becomes the crucial 
resource that allows to generate the transformation. 

For weak fields, the number of pathways connecting two states in the subspace is linearly 
related to the number of auxiliary states on the excited surface. Practically the bandwidth of 
the pulse determines this number. This means that the bandwidth in a weak field implemen- 
tation of a unitary transformations has to increase exponentially when the number of levels 
N increases. The picture is completely altered when the intensity is allowed to increase. 
Rabi cycling increases the number of interference pathways exponentially. The number of 
Rabi cycles can be estimated from the integrated intensity Jn a u ~ T/2-n Cf. Eq. (|5T|). which 
leads to an estimation of the number of interference pathways 0(M jRabi ) ~ 0(M x l 2lT ). This 
estimation is consistent with the results of Sec. II VI where only a moderate increase in X was 
observed when the number of qubits in the transformation increased. The estimated number 
of Rabi cycles changed from Jrcu ~ 6 for Q = 1 to jRabi ~ 8 for Q = 2 to Jrom ~ 12 for 
Q = 3. This means that the increase in resources of implementing a unitary transformation 
with Q qubits in a molecular environment will scale with a low power of TAE where AE is 
the pulse energy. 
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In summary, 

• A unified approach for obtaining the field that implements a unitary transformation 
has enabled the assessment of various formulations. In addition, a new algorithm 
based on the square modulus of r was developed. This scheme was found to have 
superior convergence properties with respect to the number of iterations. 

• A unitary transformation could be implemented in a molecular environment in a time 
scale of picosecond with reasonable bandwidth and intensity. For intense filed condi- 
tions the physical resources scale moderately with the number of qubits in the trans- 
formation. 

• The inversion problem of finding the field that induces a unitary transformation seem 
to be a hard numerical problem scaling unfavorably with the number of levels in the 
transformation. 
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APPENDIX A: THE VARIATIONAL METHOD 



An alternative to the Krotov method of optimization is the variational method |ll \21\. 



i. 



This methoc 
formulation 



has been used previously in the simultaneous N state-to-state transitions 
12^ and for the evolution operator formulation using the functional F re jll|. 
In the last case the variational method was generalized in terms of the evolution equation 
for the unitary transformation. Unlike the Krotov method the variational method does 
not offer a direct algorithm to minimize F sm . For simplicity only the optimization of the 
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functional F ss is discussed. The variational method is based on the functional ll 

N r T 



K({^ u ,^ fl },Ae) = J2 \(MT)\0 



A, 



o 



o s(t) 



\Ae\ 2 dt 



-2 Re 



N 



1=1 



J2(MT)\d\l) I (i> fl (t)\( Tt + -H(e + Ae))\Mt)) 



d 



i . 
K 



(Al) 



with the additional condition \ipu(t = 0)) = \l). The set of states {|/)} and the target 
unitary transformation O were introduced in section II. {\tpu(t)} denotes the initial states 
driven by the field to the final states 0|Z). The terms \^>fi{t)) are interpreted as Lagrange 
multipliers used as a constraint to impose the Schrodinger equation. The two first terms are 
equivalent to the functional (jl4J) of the Krotov method. The parameter Ao is now interpreted 
as a Lagrange multiplier. The functional (jAl|) differs from the common formulation of OCT 
in the form of the field term e + Ae. e is a reference field and Ae must be interpreted 
as the correction used to converge to the optimal field that implements the target unitary 
transformation. Setting e = and interpreting Ae as the field the common form is re- 
attained. 

By applying the calculus of variations, requiring 5K = 0, with respect to each element of 
the set {ipu(t)}, the evolution equations are reconstructed 



d 



Mt)) 



:H(e + Ae)|V«(*)> 



(A2) 



dt ft 

with the condition \ipu(t = 0)) = |/) and formal solution \ipu(t)) = U(t, 0; e + Ae)|Z). The 
variations with respect to the set {i/ifi(t)} gives 



d 



liM*)> 



:H(e + Ae)|^(*)>, 



(A3) 



dt h 

with the condition \ipfi(t = T)) = 0\l). The formal solution is \il>fi(t)) = U(£, T; e + Ae)0|/). 
Finally, variations with respect to Ae lead to the correction to the field 

N 



Ae(i) 



An H 



Im 



k (/|6 t U t (t,T;e + Ae)/2U(t,0;e + Ae)|/) 



i=i 



with 



(/|U'(T,0;e + Ae)0|/) . 



(A4) 



(A5) 



The correction to the field (|A4|) is the starting point of the iterative algorithms to find the 
optimal field. In such a case the correction to the field is implicit in the backwards and 



30 



brwards propagation of the states in Ae. Several iterative methods have been proposed 



22j |. In the simplest approach, a guess field is used to evaluate Ae, that will be used 



to obtain the input field in the next iteration. Usually it does not converge. An alternative 



procedure 



22 1 is to evaluate in Eq. (jA4j) using the field in the previous iteration and 



then to simultaneously obtain the correction to the field and evaluate U with the new field. 
This iterative algorithm is identical to the one obtained from the Krotov method in Sec. 
TTT1 A study comparing different iterative algorithms based on the Krotov and variational 
methods for the problem of state-to-state optimization is described in Ref. \2(\ . 
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